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Abstract 

When a liquid freezes, a change in the local atomic structure marks the transition to the crystal. 
When a liquid is cooled to form a glass, however, no noticeable structural change marks the glass 
transition. Indeed, characteristic features of glassy dynamics that appear below an onset tem¬ 
perature, To, m are qualitatively captured by mean field theory me], which assumes uniform 
local structure at all temperatures. Even studies of more realistic systems have found only weak 
correlations between structure and dynamics [MI]. This raises the question: is structure impor¬ 
tant to glassy dynamics in three dimensions? Here, we answer this question affirmatively by using 
machine learning methods to identify a new field, “softness?” which characterizes local structure 
and is strongly correlated with rearrangement dynamics. We find that the onset of glassy dynam¬ 
ics at To is marked by the onset of correlations between softness (i.e. structure) and dynamics. 
Moreover, we use softness to construct a simple model of slow glassy relaxation that is in excellent 
agreement with our simulation results, showing that a theory of the evolution of softness in time 
would constitute a theory of glassy dynamics. 


t Equal contribution. 
* schsani@sas.upenn.edu 
I ajliu@sas.upenn.edu 


1 



To look for correlations between structure and dynamics, one typically tries to find a 
quantity that encapsulates the important physics, such as free volume, bond orientational 
order, locally preferred structure, etc. In contrast to this approach, we use a machine 
learning method designed to hud a structural quantity that is strongly correlated with 
dynamics. Earlier, we applied this approach to the simpler problem of classifying particles 
as being “soft” if they are likely to rearrange or “hard” otherwise IE]. We describe a 
particle’s local structural environment with M = 166 “structure functions” [13] that respect 
the overall isotropic symmetry of the system and include radial density and bond angle 
information. We then dehne an M-dimensional space, with an orthogonal axis for 
each structure function. The local structural environment of a particle i is thus encoded 
as a point in M-dimensional space. We assemble a “training set” from molecular dynamics 
simulations consisting of equal numbers of “soft” particles that are about to rearrange and 
“hard” particles have not rearranged in a time Tq, preceding their structural characterization, 
and hnd the best hyperplane separating the two groups using the support vector machines 
(SVM) method [HI [H]. Finally, we dehne the softness St, of particle i as the shortest 
distance between its position in and the hyperplane, where S'* > 0 if i lies on the soft 
side of the hyperplane and Si < 0 otherwise. 

We study a 10,000-particle 80:20 bidisperse Kob-Andersen Lennard-Jones glass [TH] in 
d = 3 at different densities p and temperatures T above its dynamical glass transition 
temperature. All results here are for particles of species A only. However, the results are 
qualitatively the same for particles of both species. At each density we select a training 
set of 6, 000 particles, taken from a molecular dynamics trajectory at the lowest T studied, 
to construct a hyperplane in We then use this hyperplane to calculate Si(t) for each 
particle i at each time t during an interval of 30, OOOr at each p and T. 

We can deduce the most important structural features contributing to softness either 
by training on fewer structure functions or by examining the projection of the hyperplane 
normal onto each orthogonal structure function axis. Both analyses yield a consistent picture 
(see supplementary information): the most important features are the density of neighbors 
at the hrst peaks of the radial distribution functions gAAij") and pABij)] these two features 
alone give 77% prediction accuracy for rearrangements. Particles with more neighbors at 
the hrst peaks of g{r) have a lower softness, and are thus more stable. These results are 
reminiscent of the cage picture, in which an increase of population in the hrst-neighbor shell 
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FIG. 1. The characteristics of the softness field, a, A snapshot of the system at T = 0.47 and 
p = 1.20 with particles colored according to their softness from red (soft) to blue (hard), b, The 
distribution of softness of all particles in the system (black) and of those particles that are about 
rearrange (red). 90% of the particles that are about to rearrange have S' > 0 (shaded region). 
None of the data included in this plot were in the training set. 


suppresses rearrangements, or the free-volume picture, in which particles whose surroundings 
are closely-packed are more stable than those with more loosely-packed neighborhoods im. 
Overall, soft particles typically have a structure that is more similar to a higher-temperature 
liquid, where there are more rearrangements, while hard particles whose structure appears 
closer to a lower-temperature liquid [IS]. 

Fig. 1 (a) is a snapshot with particles colored according to their softness. Evidently, S 
has strong spatial correlations. Fig. 1 (b) shows the distribution of softness, P{S), and the 
distribution of softness for particles just before they go through a rearrangement, P{S\R). 
We see that 90% of the particles that undergo rearrangements have S' > 0. We have also 
tested other sets of structure functions (see supplemental information) and found nearly 
identical accuracy. Softness is therefore a highly accurate predictor of rearrangements that 
is reasonably robust to the set of structure functions chosen. 

We next show that the probability that particles rearrange is a function of their softness. 
This probability is calculated as the fraction of particles of a softness, S, that are rearranging 
at a given time, Pr{S). We plot Pr{S) in Fig. 2 (a) in solid lines at temperatures ranging 
from T = 0.47 (blue) to T = 0.58 (red). At each T we see that Pr{S) is a strong function 
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of softness, increasing by several orders of magnitnde, especially at the lower temperatnres, 
in the range S = —3 to S' = 3. A similar, albeit more modest, relationship was seen 
in [in]. When Pr{S) is plotted as a fnnction of 1/T for several values of softness. Fig. 2 
(b), the probability that a particle of softness S will rearrange has Arrhenius behavior, 
Pr{S) = Po{S) exp{—AE{S)/T) where Po{S) and AE{S) depend on S. Conhrming this 
observation, Pr{S)/Po{S) collapses over many orders of magnitude for all temperatures 
when plotted against AE{S)/T, as shown in the inset of Fig. 2(b). 

An Arrhenius form emerges when a kinetic process depends on a single energy scale 
AE{S). In Fig. 2 (c) we plot AE{S) and F(S') = lnPo(*S') vs. S. Both terms depend nearly 
linearly on S: AE = cq — eiS and S = Sq — SiS* where all four coefficients are positive and 
independent of T. Our results are consistent with the interpretation that at low tempera¬ 
tures, harder regions of the glassy liquid with higher energy barriers are frozen out while 
softer regions are not, leading to heterogeneous dynamics. These heterogeneities smooth 
out with increasing temperature, and vanish altogether once Pr{S) no longer depends on 
softness. This occurs at the temperature Tq where the softness dependence of S exactly 
cancels that of AE/Tq and so Tq = Ci/Si. This result can also be seen visually in Fig. 2 
(b) where the different Arrhenius predictions for Pr{S) all intersect at a single temperature. 
To, where the probability of rearrangement will be independent of softness. In Fig. 2 (d) 
we compare our prediction for Tq to the onset temperature of glassy dynamics measured by 
Keys et al. [20], T™, at different densities. The excellent agreement between the predicted 
To and the measured values implies that the onset of glassy dynamics at T = To coincides 
with the onset of correlations between structure (softness) and dynamics. 

We explore next the relationship between softness and the nonexponential decay of the 
overlap function 

(lit) = 

i 

where N is the number of particles in the system, r* is the position of particle i, and 0 is 
the Heaviside function. We take a = 0.5 [21]. In Fig. 3 (a)-(b) we plot the overlap function 
for different temperatures at p = 1.20. Our aim is to understand the form of the decay of 
q{t) from the behavior of the rearrangement probability, Pr{S). To begin, we dehne the 
contribution to the overlap from particles whose softness was initially S' at t = 0, q{S,t). 
The total overlap is q{t) = f dSq{S,t)P{S). Because qiS,t) is the fraction of particles with 
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FIG. 2. The relationship between softness and dynamics, a, The probability that par¬ 
ticles rearrange as a function of their softness, Pr{S), for temperatures T =0.47, 0.53, and 0.58 
plotted in blue to red. Solid lines are measurements from molecular dynamics trajectories (solid 
lines). Dashed lines present the probability computed using the Arrhenius form for Pr{S) (dashed 
lines). Points represent the probabilities calculated from the zero-time derivative of the overlap, 
—dq{S,t)/dt at T = 0.47 and T = 0.58. b, Pr{S) as a function of 1/T for 5 different softness 
values from S ~ — 3 (blue) to S ~ 3 (red). The inset shows the collapse of these probabilities when 
Pr/Pq is plotted against AE/T. c, AE and S, where Pr{S) = exp(S — AE/T), vs. softness S. 
d, predicted onset temperature Tq vs. T™, onset temperature measured by Keys, et al. [20j, for 
densities p = 1.15,1.20,1.25,1.30. The straight line corresponds to Tq = T™. 


initial softness S that have not rearranged after a time t, we expect = —CaPR{S) 

(see Supplementary Information for details) where Ca is the fraction of rearrangements that 
displace particles by more than a. This is indeed the case, as is evident from the data in 
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FIG. 3. Overlap calculated from softness a, Solids lines are the measured overlap function, 
for temperatures T = 0.45, 0.47, 0.53, 0.58, 0.63, and 0.70, from blue to red, respectively. The 
dashed lines in the insets show predictions assuming each Arrhenius process is independent of one 
another, b. The solid lines in the insets are the same as in a. Dashed lines are predictions for the 
overlap function from Pr{S) including changes in the softness field induced by spatial correlation 
between rearranging particles. 

Fig. 2 (a), when |f=o (points) is overlaid with Pr{S) (solid lines). 

If we now assume that each particle rearranges with probability Pr{S) as an independent 
Arrhenius process according to Fig. 2, then we can predict the decay of q{S, t) using a 
simple discrete model; it can be written in terms of the probability that a particle of 
softness S does not rearrange for t — 1 timesteps before dually rearranging at time t, (1 — 
PR{S)y~^P r{S) . The resulting prediction for q(t) (dashed) is shown in Fig. 3 (a) for several 
different temperatures. While the prediction is not poor, its accuracy decreases at longer 
times, particularly at lower temperatures. 

We now show that the discrepancy between our naive theory and the decay of q{t) pri¬ 
marily results from a crucial neglected feature: even if a given particle does not rearrange, 
its local structural environment -and therefore its softness- can be altered if nearby particles 
rearrange. This physics is reminiscent of facilitation. 

To take this facilitation into account, we calculate the “softness propagator”, G{S, So,t), 
that is, the distribution of softness at time t for particles that start with a softness Sq at 
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FIG. 4. Time evolution of softness a, The stochastic evolution of softness in time as seen in 
through the evolution of the Gaussian approximation to the distribution of softness, b, The time 
evolution of the softness distribution for a collection of particles with initial softness ~ — 3 from 
t = 0 (blue) to t = lOOOr (pink). Points are the measured histogram values, and the dashed lines 
are Gaussian approximations to the distribution, c, The time evolution of the average softness for 
particles that start from several softness values ranging from ~ — 3 (blue) to 5o ~ 3 (red). 


t = 0 and move less than a distance a after a time t {i.e. that do not rearrange in a time 
t). Fig. 4 (a) shows a Ganssian approximation to G{S, Sq = —3,t). We see that G{S, So,t) 
is sharply peaked aronnd Sq at small t but widens and shifts with increasing t reminiscent 
of directed diffusion. Fig. 4 (b) shows G{S,So = —3,t) at several different times, where 
points are measured probabilities and dashed lines are their Gaussian approximations. In 
Fig. 4 (c) we plot the mean softness evaluated as {S{t))so = f dSSG{S, So,t) for several 
different values of Sq. For each Sq the average softness of particles evolves towards the mean 
of the equilibrium softness distribution over a time period of approximately r^. The softness 
propagator is evaluated only for particles that have not rearranged, so Fig. 4 shows that 
rearrangements of nearby particles affect a particle’s softness signihcantly. 

Our hrst naive prediction based on the assumption that particles rearrange independently 
corresponds to G{S,SQ,t) = 6{S — Sq). We rehne our theory by using the actual softness 
propagator in connecting the probability of rearranging, Pr(S'), with the overlap q{S,t) 
(see Supplementary Information). For ease of calculation, we approximate G{S,SQ,t) as a 
Gaussian distribution in S and calculate its mean and variance as functions of Sq and t 
from simulated data. The resulting prediction for the overlap is shown in Fig. 3 (b). The 
agreement with the actual q{t) is excellent, suggesting that an understanding of the time 
evolution of the softness held, or equivalently of the softness propagator, would suffice to 
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understand the non-exponential decay of the overlap function. 

Our results show that there is structure hidden in the disorder of glassy liquids. This 
structure can be quantihed by softness, which controls glassy dynamics at temperatures 
below Tq. According to our analysis, simple Arrhenius relaxation for each softness, coupled 
with the time evolution of softness, leads to the observed slow, non-exponential relaxation 
dynamics of glassy liquids below Tq. Thus, our results suggest that the challenge of under¬ 
standing glass transition dynamics can be reframed as the challenge of understanding the 
evolution of softness. 
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discussions. This work was supported by the U.S. Department of Energy, Office of Basic 
Energy Sciences, Division of Materials Sciences and Engineering under Award DE-FG02- 
05ER46199 (S.S.S. and A.J.L.), and Harvard lACS Student Scholarship (E.D.C.). This work 
was partially supported by a grant from the Simons Foundation (7^305547 to A.J.L.). 


I. METHODS 

A. System information. 

We study a 10,000-particle Kob-Andersen model, a 80:20 binary LJ mixture [16] with 
parameters: aAA = 1-0, (Jab = 0.8, asB = 0.88, e^A = 1-0, (ab = 1-5, (bb = 0.5. Time 
is measured in units of r = \/(aa!(j\a Boltzmann constant is ks = 1. We cut 

off the LJ potential at 2.baAA and smooth the potential so that force varies continuously. 
This mixture has been characterized extensively. In particular, we compare our predictions 
to the measurements of the onset temperature in Keys et al. [20]. Simulations were done 
using LAMMPS [ 22 ] in an NVT ensemble with a Nose-Hoover thermostat and a timestep 
of 0.0025r. We output states every r and quench them to their nearest inherent structure 
using a combination of conjugate gradient and FIRE algorithms. Throughout this study 
we use inherent structure positions. However, qualitatively similar results can be obtained 
using time averaged positions. We study this system over the temperatures and number 
densities listed in Table III 



TABLE I. Number densities and temperatures studied. Each column contains the temperatures 
studied for a given number density p. 


p 

1.15 

1.20 

1.25 

1.30 

T 

0.37 

0.47 

0.58 

0.70 


0.42 

0.51 

0.61 

0.75 


0.45 

0.53 

0.69 

0.84 


0.52 

0.56 

0.76 

0.92 



0.58 

0.97 

1.12 



0.70 

0.97 



B. Identifying rearrangements. 

We adapt a method first proposed by Candelier et a/. [23l El]. A timescale tR = lOr 
is chosen to be commensurate with the amount of time the system takes to complete a 
rearrangement. Then two time intervals are dehned as A = \t — tR/2,t] and B = \t,t + tR/2]. 
An indicator function can then be written as, 

Thop(t) = \/ {{fi - {ri)B?)A{{fi - {ri)A?)B (1) 

where ()a and ()b are averages over the intervals A and B respectively, phop is large when 
the mean position of a particle changes appreciably. Otherwise, it is similar in magnitude 
to the variance in particle positions due to noise from the inherent structure calculation. 

To hnd rearrangements we restrict our attention to events in which p^op exceeds a thresh¬ 
old of 0.05, that is large compared to the scale of fluctuations in particle positions but 
small compared to the typical value of phop during a rearrangement. As discussed in the 
supplementary material, we define rearrangements to be those events with pj^^p > Pc = 0.2. 
Changing this cutoff affects the results only quantitatively and manifests itself primarily as 
a shift in the energy scale, that is approximately logarithmic in the cutoff. This agrees 
with the observations of Keys et al. [201 who saw a similar logarithmic shift in the energy 
scale governing rearrangements with the size of the rearrangements. 

Note that rearrangements defined using phop result in particle displacements that follow 
a distribution that depends on the cutoff pc used. This pc dependence needs to be addressed 
when comparing the probability of rearrangement to the overlap function and its derivative. 
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which are defined in terms of a length scale a. To do this, we mnltiply Pr by a temperature- 
independent constant Ca, namely the fraction of rearrangements that displace particles by 
more than a. 


C. Computing softness. 

We have made two improvements that greatly increased the prediction accuracy for re¬ 
arrangements compared to Ref. [12]. First, we identified rearrangements more carefully, as 
detailed above. Second, we defined our training sets more carefully. Each training set con¬ 
tains 6000 particles that rearrange in the next time step, each labeled with r, = 1, as well 
as 6000 particles that have not rearranged for a time Tq, before the structure was calculated, 
each labeled with rj = 0. These particles were chosen randomly from the set of all particles 
satisfying these conditions from MD simulations at a low temperature. Then, a training set 
of N particles can be written as {(Fi,ri) tn)}, where Fi = are the 

the M structure functions that describe the local neighborhood of particle i [12]. We then 
use an SVM to find the hyperplane w ■ F — b = 0 that separates the points with r* = 1 
from those with = 0. This hyperplane is used on the rest of the data to reach the results 
reported. 

The SVM is trained, that is, the hyperplane is constructed, on the binary variable r using 
the LIBSVM package [H]. It is not possible to find a hyperplane that perfectly separates 
the two different classes. We use a penalty parameter C and find the optimal hyperplane 
equation by minimizing 

1 ^ 

-w'^ -w + C^^i, ( 2 ) 

i=l 

with the constraint yt ■ {w'^ ■ Fi + h') > 1 — and > 0. The C parameter was chosen 
through cross-validation [T2]. The hyperplane obtained from this training can be used to 
classify a new particle neighborhood, Fn, as soft or hard. Fn is soft if m ■ Fk — 6 > 0, and 
hard otherwise. The continuous variable softness is defined hj Sn = w ■ Fn — b. Training 
a neural network to classify soft and hard particles, and using the output from the hidden 
layer of the neural network as softness, yields similar results. Here we use only the SVM 
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approach. 
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